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ABSTRACT 

The goal of this paper is to derive signatures of synchrotron radiation from state-of-the-art simulation 
models of collimating relativistic magnetohydrodynamic (MHD) jets featuring a large-scale helical 
magnetic field. We perform axisymmetric special relativistic MHD simulations of the jet acceleration 
region using the PLUTO code. The computational domain extends from the slow magnetosonic 
launching surface of the disk up to 6000 2 Schwarzschild radii allowing to reach highly relativistic 
Lorentz factors. The Poynting dominated disk wind develops into a jet with Lorentz factors of T ~ 8 
and is collimated to 1°. In addition to the disk jet, we evolve a thermally driven spine jet, emanating 
from a hypothetical black hole corona. Solving the linearly polarized synchrotron radiation transport 
within the jet, we derive VLBI radio and (sub-) mm diagnostics such as core shift, polarization 
structure, intensity maps, spectra and Faraday rotation measure (RM), directly from the Stokes 
parameters. We also investigate depolarization and the detectability of a A 2 -law RM depending on 
beam resolution and observing frequency. We find non-monotonic intrinsic RM profiles which could 
be detected at a resolution of 100 Schwarzschild radii. In our collimating jet geometry, the strict 
bi-modality in polarization direction (as predicted by Pariev et al.) can be circumvented. Due to 
relativistic aberration, asymmetries in the polarization vectors across the jet can hint to the spin- 
direction of the central engine. 

Subject headings: galaxies: active - galaxies: jets - ISM: jets and outflows - magnetohydrodynamics 
(MHD) - radiation mechanisms: non-thermal - relativistic processes 



1. INTRODUCTION 

Relativistic jets are launched from accretion disks 
around compact objects and accele rated and collimated 
by m agnetohydrodynamic forces (Blandford k Payne 
1982[ ). The disk jet flow may be accompanied by a 
highly relativistic central spine jet resulting from elec- 
trodynamic effects within t he black hole magnetosphere 
( [Blandford k Znajek|[l977| . 

The formation of relativistic MHD jets has been in- 
vestigated by a number of authors since the seminal pa- 
per of Blandford k Znajek. The early attempts were to 
look for stationary state MHD solutions for the asymp 
totic jet structure (Chiueh et al. 1991 Appl k Camen 



zind 



of co 



1993 Fendt 



1997b), or the jet formation domain 



limation and acceleration in case of self-similarity 



(Contopoulos 1994 1995), or taking into acc ount the 



2.5D force- balance in case o f spec ial relativity (Camen- 



f ali|[l98 7; Fendt k Mcmola 2001 }, or genera l relativity 
todt|l997a| |Ghosfi||2000[ |Fendt k Greiner||2001| ). 
Time-dependent simulations of jet formation from the 
disk surface have first b een investigated in the non- 
relativistic appr oximation ( Ustyugova et al.|1995 Ouyed 
k Pudritz 1997). General relativistic MHD simulations 
ot accretion disks launching outflows (INishikawa et al. 
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2005; McKinne^ 
McKiimey k 



ley fc Nara yan 2007] jBeckwith et al.pOOSl 
Blandford 2009J), indicate that highly rela- 
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tivistic jets may not be launched by the disk itself, but 
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from the black hole magnetosphere. These jets - so-called 
funnel flows - may reach Lorentz factors up to 10 how- 
ever, the mass loading is put in by hand (maybe corre- 
sponding to a pair-creation process) and is not relat ed to 
the accretion process itself. |Komissarov et al. ( 2007 ) pro- 
posed special relativistic MHD simulations of AGN jets 
over a large spatial scale and find asymptotic Lorentz 
factors of about 10. These jets were pressure-confined 
by the outer boundary condition. 

In a recent paper we have treated another s etup of rel- 
ativistic jet formation (Porth & Fendt 2010). We have 
applied special relativistic MHD to launch outflows from 
a fixed-in-time hot surface of an accretion disk which is 
rotating in centrifugal equilibrium with the central com- 
pact object. This allowed us to pass through the slow- 
magnetosonic point, thus obtaining consistent mass and 
Poynting fluxes, and to investigate the subsequent accel- 
eration and collimation degree of a variety of jets. Es- 
sentially, we find these relativistic MHD jets to be truly 
self-collimating - similarly to their non-relativistic coun- 
terparts. In this paper, we will prescribe certain mass 
fluxes and explore the influence of the poloidal current 
distribution on the jet formation process. In contrast 
to the previous work, the jet energy and its partitioning 
is now a true parameter of the models, allowing us to 
investigate the regime of highly relativistic flows. More 
importantly, observational signatures of the jet formation 
models are discussed. 

The existence of an ordered large-scale /zG- mG mag- 
netic field in extragalactic jets is well established by the 
detection of radio synchrotron emission, however the ex- 
act geometry of the field structure cannot be derived. 
The recent observational literature strongly suggests he- 
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lical magnetic fields (e.g. O'Sullivan & Gabuzda 
a scenario which is consistent with theoretica. 



(2009)), 
models 



and numerical simulations of MHD jet formation (see 
e.g. dBlandford fc Payne||1982[ |Ouyed fe Pudritz||1997| 
' Krasnopolsky et al.||1999| |Fendt fe Cemeljic|2002| (Fendt 
2006| |Porth fe Fendt | |2010p ). The wavelength dependent 



have taken a step into this direction by providing opti 
cally thin synchrotron a nd polarization maps from self- 



similar MHD solutions. Broderick & McKinney (20101 



rotation ot the polarization plane known as Faraday ro- 
tation provides a valuable diagnostic of magnetic field 
structure in astrophysical jets. Consistent detections of 
A% oc A 2 are foun d in resolved jets as w ell as in unre- 
solved radio cores ( |Zavala fe Tay lor 2003). Helical mag- 
netic fields are generally perceived to promote transver- 
sal Faraday rotation measure (RM) gradients owing to 
the toroidal field component. Observationally, such gra- 



dients were first detect ed by Asada et all (|2002|) and 

7273. 



Zavala fc Taylor| ( |2005[ ) in the jet of 3C 273. The RMs 
are generally found to follow a monotonic profile across 



jets (Gabuzda et al. 2004 O'Sullivan & Gabuzda 2009 
201U|) 



Croke et al 

While the magnetic field structure of pc and Kpc-scale 
jets can in principle derived from radio observations, not 
much is known about the field structure within the jet 
forming region very close to black hole - mainly because 
of two reasons. Firstly, this very core region cannot be 
spatially resolved and, thus, cannot readily be compared 
with expected RM profiles of helical jets. This is some- 
what unfortunate as close to the jet origin we expect 
the magnetic field helix to be well preserved and not 
much affected by environmental effects. Secondly, the 
observed rotation measures are so high that it is impos- 
sible to draw firm conclusions about the intrinsic field 
geometry from the polarization vector with mere radio 
observations. 

Only very few cases exist where this core of jet forma- 
tion could be resolved observationally. Among them is 
the close-by galaxy M87, where the V1BI/V1BA reso- 
lution of ~ 0.1 mas is sufficient to resolve about 0.01 pc 
within the central region, and allows to trace the jet ori- 
gin down to ~ 20 Schwarzschild r adii rs when the recent 
mass estimate of 6.6 x 10 9 M Q by Gebhardt et al (2011 1 
is adop ted. This pinpoints the launching ar ea within 
~ 30r s ( |Junor et allfll999| jKovalev et al.||2007[ ). The ra- 
dio maps clearly show limb brightening and indicate an 
initial jet opening angle of about 60°. 

Despite a vast amount of observational data spanning 
over a huge frequency scale, and also time series of these 
multifrequency observations, rather little is known about 
the dynamical status of relativistic jets. There are no di- 
rect unambiguous observational tracers of jet velocity or 
density as only (if at all) the pattern speed of radio knots 
is detected. Kinematic modelling of knot ejection sug- 
gests pc-scale lorentz factors of typically T ~ 10, while 
Kpc-jet velocities are believed to be definitely lower and 
of the order of 0.1c. Kinematic modeling of jet prop- 
agation has been combined with synchrotron emission 
models of nuclear flares, resulting in near perfect fits of 
the observed, ti me-dependent radio p attern of jet sources 
such as 3C279 ( jlindfors et aT|2006|) 

However, what was missing until recently is a con- 
sistent combination of dynamical models with radiation 
models of synchrotron emission resulting in theoretical 
radiatio n maps which can then be c ompared with obser- 
vations. Zakamska et al. (2008) and Gracia et al. (20091 



presented synchrotron ray-castings from 3D general rela- 
tivistic jet formation simulations that also include the 
evolution of a turbulen t accr etion disk performed by 
McKinney & Blandford (2009). In their approach, the 
MHD solution is extrapolated by means of an essentially 
self-similar scheme in order to reach distances up to lOpc. 
Their study focussed on the rotation measure provided 
by Faraday rotation in the disk wind external to the emit- 
ting region in the Blandford-Znajek jet. 

In comparison, using axisymmetric large scale simu- 
lations, we do not rely on an extrapolation within the 
AGN core (up to 0.3pc) and treat the Faraday rotation 
also internal to the emitting region in the fast jet that 
gradually transforms into a sub-relativistic disk wind. 
The observational signatures obtained in our work are 
derived entirely from the (beam convolved) Stokes pa- 
rameters which allows us to also investigate the breaking 
of the A 2 law due to opacity effects. 
As the above studies, we also rely on post-hoc prescrip- 
tions for the relativistic particle content. Towards a more 
consistent modeli ng of the non-thermal particles, |Mim-| 
ica & Aloy (2010) have presented a method to follow the 
spectral evolution of a seed particle distribution due to 
synchrotron losses within a propagating relativistic hy- 
drodynamic jet. The question of particle acceleration 
and cooling is in fact essential to close the loop for a 
fully self-consistent treatment of jet dynamics, jet inter- 
nal heating, and jet radiation. 

To derive signatures of synchrotron emission from rela- 
tivistic MHD jet formation is the main goal of the present 
paper. We apply the numerically derived dynamical vari- 
ables such as velocity, density, temperature, and the mag- 
netic field strength & configuration to calculate the syn- 
chrotron emission from these jets, taking into account 
proper beaming and boosting effects for different incli- 
nation angles. In particular, we apply the relativistically 
correct polarized radiation transfer along the line of sight 
throughout the jet. 

2. THE RELATIVISTIC MHD JET 

We perform axisymmetric jet acceleration simul ations 
with the P1UTO 3.01 code flMignone et~aT||2007[ ) solv- 
ing the spe cial relativistic magne tohydrodynamic equa- 
tions. As in|Porth & Fendt (20101, the simulations are of 



the disk-as-boundary type, however in the present study 
the jet starts out with the slow magnetosonic velocity 
opening the freedom to assign all energy channels as 
input parameters. The energy of the jet base is dom- 
inated by Poynting flux, driving a large-scale poloidal 
current circuit. This current distribution is prescribed 
as a boundary condition in the toroidal magnetic field 
component. We investigate three cases for B${r) oc r~ s 
with s G {1, 1.25, 1.5}, where r is the cylindrical radius, 
resulting in asymptotic jets which are either in a current- 
carrying (one of them) or current-free (two others) con- 
figuration. Injected initially with slow magnetosonic ve- 
locity, the jet material accelerates through Alfven and 
fast magnetosonic surfaces within the simulated domain. 

2.1. Numerical grid setup 



Jet base synchrotron radiation 
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To eradicate any artificial collimation effect from the 
outflow boundaries, we decided to move the domain 
boundaries to such large distance that they are out 
of causal contact with the region of interest. An in- 
ner equidistant grid of 20 cells extends to the scale ra- 
dius 7"i corresponding to the inner radius of the accre- 
tion disk. Beyond T\ the grid is linearly stretched to 
(r, z) = (500,500) with a scaling factor of L0013, and 
for (r, z) > (500, 500) with a larger scaling factor of 
f.0047. So far, we apply a square box of 2000 2 scale 
radii corresponding to 2555 2 grid cells. Magnetic fields 
are advanced on a sta ggered grid using the m ethod of 
constrained transport (Balsara & Spicer 19991 supplied 
by PLUTO. 

For the calculation of the radiation maps, we consider 
only a subdomain of 1417 x 2356 grid cells corresponding 
to 200 x 1000 scale radii. 

The first possible contamination by boundary effects 
takes place when the bow shock reaches the upper out- 
flow boundary after t cl < z cn d/c = 2000. A typical sim- 
ulation is terminated after maximally i c2 = 3000 time 
units when a signal traveling at the speed of light could 
have returned from the z en d boundary to the subdomain. 
This conservative treatment ensures that no spurious 
boundary effects can occur. 

In order to test this we have re-run one of our sim- 
ulations on a five times larger computational domain 
with comparable resolution. This simulation provided 
the same result as the lower grid size simulation and ac- 
quired a nearly stationary state when terminated. Thus, 
the region of interest is not at all affected by any outflow 
boundary effects on account of an increased computa- 
tional overhead. 

2.1.1. Inflow boundary conditions 

In a well posed MHD boundary, the number of out- 
going waves (i.e. seven minus the downstream critical 
points) must equal the number of boundary constraints 
provided. Thus, in addition to the B^(r) and v p (r) pro- 
files, we choose to prescribe the thermal pressure p{r) 
and density p(r) as boundary conditions for the jet in- 
jection (the jet inlet). The fifth condition sets v P ||-B p 
and thus constrains the toroidal electric field E${r) = 0. 
This is a necessary condition for a stationary state to be 
reached by the axisymmetric simulation. 

The remaining primitive MHD variables and B r are 
extrapolated linearly from the computational domain, 
while the component B z (which determines the magnetic 
flux) follows from the solenoidal condition. 

Specifically, the fixed profiles read 



p{R)=Pi [(1-9)R + 6R- 1 - 5 ] 
p(R) = Pl [(1 - 0)(i - px lni?) + 9R- 2 - 5 ] 
B^r) = B^ [(1 - 6)r + 9r~ s ] 
v p {r) = v sm {r) 



where 



r < 1 
r > 1 



(1) 

(2) 
(3) 
(4) 

(5) 



D sm (r) of relation 
count for the varia 



ing, a subscript 1 indicates the quantity to be evaluated 
at (r,z) = (1,0). 
The slow magnetosonic velocity profile along the disk 
4| is updated every time step to ae- 
ries which are extrapolated from the 
domain. Within R < 1, an inner "black hole corona" 
with relativistic plasma temperatures is modeled. The 
physical processes responsible for the formation of an 
inner hot corona could be an accretion shock or a so- 
called CENBOL shock (CENtrifuga l pressure supported 
BOundary Layer shock) (see also |Kazanas fc E llison 
19861 l Das fc Chakrabarti|[l999| . Also |Blandford| 
has proposed a mechanism of dissipation near the ergo- 
sphere as a consequence of the Lense-Thirring effect. Due 
to the large enthalpy and decreasing Poynting flux of the 
inner heat bath, a the rmally driven outflo w is anticipated 
from this region (e.g. |Meliani et aL|2010|) We apply the 
causal equation of state introduced by |Mignone et al.| 
(20051 to smoothly join the relativistically hot central 
region to the comparatively cold disk jet. This choice 
permits a physical solution for both regions of the flow 
(see also Mignonc fc McKinney1[2007l ). 



2.1.2. Outflow boundary conditions 

At the outflow boundaries we apply power-law extrap- 
olation for p, p and the parallel magnetic field component, 
while we apply the solenoidal condition to determine the 
normal magnetic field vector B z respectively B r . For the 
velocities and the toroidal field, zero-gradient conditions 
are applied. This choice is particularly suited to preserve 
the initial condition that is well approximated by power 
laws. A more s ophisticated treatment s uch as force-free 
(introduced by SRoma nova et aL] 1 1997 1 or zero-current 
boundaries as in |Forth fc Fend" ( |2010 ) is rendered un- 
necessary by the increased computational grid as detailed 
before. 

2.2. Initial conditions 

As initial setup we prescribe a non-rotating hydro- 
static corona threaded by a force-free magnetic field. 
For the initial poloidal field distribution we adop t B r — 



l/r-z/r ( r 2 + z 2 ) 1/2 an d B z = (r 2 + z 2 ) 1/z ( |Ouyed 



1/2 



is the step function and r and R denote the cylindrical 
and spherical radius, respectively. Here and in the follow- 



fc Pudritz| ( [19971 ), see also l Porth fc Fei^dtj pHPf ). The 
hydrostatic corona is balanced by a point-mass gravity in 
a Newtonian approximation. In order to avoid singulari- 
ties in the density or pressure distribution, equations (JlJ), 
([2j) , we slightly offset the computational domain from the 
origin by (r , z ) = (0, 1/3). 

2.3. Parametrization 

In the present setup we focus on effects of the poloidal 
current distribution and parametrize accordingly. Thus, 
we set the Kepler speed at r — 1 to vk — 0.5c for conve- 
nience, yielding a sound speed c Sj i = (7Pi/pi) 1 ^ 2 — 0.4c 
from the hydrostatic condition. To further minimize 
the number of free parameters, we tie the toroidal field 
strength given by B^^ to the poloidal field strength via 
B<j, t i — 0.5B p i. This Ansatz is consistent with the sub- 
Alfveninc nature of the flow, since at the Alfven point 
B^i ~ B Pi \ (e.g. Krasnopolsky et al.|1999 ) is valid. The 
two remaining parameters are the poloidal magnetic field 
strength measured by the plasma-beta (if — 2px/B 2 A 
and the toroidal field profile power law index s. 
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Note that with the choice of a fixed in time toroidal 
field, the injected Poynting flux is controlled via the /3i 
parameter, implying the toroidal field being induced in 
a non-ideal MHD disk below the domain. The conserved 
total jet energy flux and its partitioning between Poynt- 
ing and kinetic energy given by the cr-parameter essen- 
tially become boundary conditions and determine the 
terminal Lorentz factor. 

2.4. Physical scaling 

In order to convert code units to physical units we 
need to define the two scaling parameters of length and 
density, while the velocity is naturally normalized to the 
speed of light c. To obtain an approximate radial scale, 
we assume that the transition between the inner corona 
to the disk-driven jet at r — 1 corresponds to the in- 
nermost stable circular orbit at 3i?s (for a Schwarzschild 
black hole). Hence the physical length scale is given by 

:m. (6) 



8.9 x 10 



14 



10 9 M G 



The physical density is then obtained by assuming a total 
jet power E± 3 in units of 10 43 erg s _1 , 



4.7 x 10 



-17 -^43 / M, 



E V 109M G 



p g cm 



(7) 



where E is the corresponding power in code units for a 



certain simulation run (of order 
magnetic field strength follows to 



10). With this, the 



B r 



7.3 x 10 2 




B Gauss 



(8) 



(9) 



and the physical time scale becomes 

Unless stated otherwise, we adopt a black hole mass 
of M. = 10 9 A/ Q and a total jet power of E43. In order 
to calculate the observable radiation fluxes, we assume 
a photometric distance of D = 100 Mpc. The angular 
scale of the Schwarzschild radius then becomes a r g = 
0.2 /xas. For the case of M87's supermassive black hole 
with M. = 6.6 x 1O 9 M and D = 16 Mpc we would 
have a r s = 8 //as yielding an increase in resolution by a 
factor of 40 compared to our fiducial scaling. 

3. JET DYNAMICS: ACCELERATION AND COLLIMATION 

The large intrinsic scales of the relativistic MHD jet 
acceleration process require substantial numerical effort 
when simulated with a dynamic al code. Codes opt i mized 
fo r such tasks were deve l oped by Komissarov et al. ( 2007 ) 
or Tchekhovskoy et al. ( 2008 ) , involving particular grid- 
extension techniques which can speed up the simulations 
of a causally de-coupled flow. In both seminal papers 
the flow acceleration could be followed substantially be- 
yond the equipartition regime to establish tight links to 
analytical calculations. However, it can be argued that 
by using a rigid, reflecting bounda ry of certain shape 
as done by Komissarov et al. ( 2007), or a forc e-free ap- 
proach as applied by Tchekhovskoy et al. ( 2008 ), the rate 



of jet collimation resulting from those simulations could 
be altered. 

Here we aim at studying MHD se//-collimation includ- 
ing inertial forces. We therefore solve the full MHD equa- 
tions omitting the outer fixed funnel around the jet and 
replace it with a stratified atmosphere which may dynam- 
ically evolve due to the interaction with the outflow. By 
placing the outer boundaries out of causal contact with 
the solution of interest, we can be certain to observe the 
intricate balance between jet self-collimation and accel- 
eration that is inaccessible otherwise. We like to em- 
phasize that the density has to be considered in the flow 
equations for two reasons - one is to take into account 
the inertial forces which are important for collimation 
and de-collimation, the other is our aim to consistently 
treat the Faraday rotation which is given by cold elec- 
trons. The latter could in principle be taken into account 
in magnetodynamic simulations by introducing electron 
tracer particles as additional degree of freedom, but is 
not immediately satisfied by applying the force-free limit 
of ultrarelativistic MHD. 

We summarize our parameter runs of different jet mod- 
els in table [TJ indicating the maximum Lorentz factor 
attained, r max and other dynamical quantities of inter- 
est on which our results discussed in the following are 
based. Although the different electric current distribu- 



Table 1 

Simulation runs 



run ID 


s 


Pi 


Tmax 


Too 


%,l 


lh 


1 


0.005 


8.5 


25 


0.21° 


2b 


1.25 


0.005 


7.9 


21 


0.16° 


3h 


1.5 


0.005 


7.9 


23 


0.17° 


1m 


1 


0.01 


5.9 


13 


0.14° 


2m 


1.25 


0.01 


5.6 


13 


0.16° 


3m 


1.5 


0.01 


5.9 


13 


0.24° 



Note. — Columns denote: Simulation ID; Radial power law 
slope of the toroidal field; Plasma beta at (r,z)=(l,0); Maximal 
Lorentz factor obtained in the Raycasting domain; Maximal at- 
tainable Lorentz factor assuming complete conversion into kinetic 
energy Too = fi. Collimation angle of the field line rooted at the 
inner disk radius evaluated at z=1000. 



tions applied in the inflow boundary condition promote 
a distinct jet dynamics as seen for example in the po- 
sition of the light cylinder shown in figures [l]to[3| the 
geometry of the field lines turns out to be quite similar. 
At a height z = 750 (corresponsing to w 2200rg), the 
fast jet component is collimated into an opening angle 
less than 1° in all our models. More obvious differences 
are found in the radial distribution of the Lorentz fac- 
tor which is peaked at the maximum of vertical current 
density j z . In the case of closed-current models, the fast 
jet component becomes narrower as the integral electric 
current levels off more steeply. We find an acceleration 
efficiency in terms of the total energy per rest mass en- 
ergy, fi of T/fi ps 80% for the axial spine and acceleration 
efficiencies varying between 20% < T/fi < 40% for the 
outer parts of the jet (see the outer field lines in the mid- 
dle panels). Since the flow has not reached equipartition 
within the considered domain, the acceleration efficien- 
cies we obtain represent only a lower limit to the total 
efficiency T^/ p. Energy conversion is depicted in the 
right panels of figures [l] to [3] showing the individual en- 



Jet base synchrotron radiation 



■5 



1000 



800 



600 



400 



200 











* - 1 5 

-20 
=g§ 








0.0 1 00 
* 0.0010 
0001 






100.0 




a. io.o 




q 


b " 1.0 
0.1 










j:e 

=*- 0.6 

?: 0.4 




8:3 




10.00 








E 1.00 

0.10 
0.01 






T~ 80 
— 60 
* 40 

*" 2 8 





- A 


' F J c Total 


\ Enthalpy 


i i 

i | Kinetic 

[ i Poynting 


--. « Gravity 







0.2 0.4 0.6 0.8 1.0 1.2 




20 40 60 80 100 



10 15 20 25 



Figure 1. Jet model lh allowing no outgoing current in the disk. Left: Current lines (dashed) shown on Lorentz factor color-contours in 
the (r,z) plane (note the extreme aspect ratio). Field lines are given in solid white and the light-cylinder is indicated by the solid black 
line. Center: Cuts through z = 750 for Lorentz-factor V, integral current I z , lab-frame density Tp, Poynting-to-kinetic energy flux ratio 
cr (solid) and total normalized energy fi (dashed), acceleration efficiency T/fi, field line collimation angle 0g and the pitch angles of the 
co- moving system (dashed), respectively the lab-frame (solid). Right: Acceleration along selected field-lines against the cylindrical 
radius r snowing thermal acceleration for a field line in the spine (footpoint rf p = 0.2, above) and magnetic acceleration in the jet (rf p = f .5, 
below). Vertical lines indicate the crossing of the Alfven (A) and fast (F) critical point as well as the light cylinder (lc). 



ergy channels normalized to the conserved rest mass en- 
ergy flux pUpC 2 along selected field lines. The terms are 

defined: -^Enthalpy = T(h-1), ^Kinetic = T, ^Gravity = V 

and the Poynting flux £p oy ntmg = — rVLBpB^ / (At: pu p c 2 ) . 
[^] In the latter relation, we have introduced Ferraro's 
iso-rotation parameter f2 defining the "angular velocity 
of the field line" given by 



I > 

B,, 



(10) 



& Fendt ( 2010[ )). Thus, the asymptotic flow acceleration 
depends solely on the decrease of 

= r 2 B p (13) 

alongthe flow line by differential fanning out of the field 
linesjj We show the evolution of the <f> function along 
selected field lines for the intermediate model 2h in figure 
[4] (left panel). In the non-asymptotic regime, <j) increases 
until the x = 1 surface, while it is decreasing for x 3> 1 
as expected. 



3.1. Poynting dominated flow 

We first consider the MHD acceleration of the disk 
component of the jet flow. Here, the bulk of the ac- 
celeration takes place in the relativistic regime beyond 
the light surface rr,(r, z) = c/f2, and can therefore be ap- 
proximated asymptotically x = r/r^ 3> 1, T 1. For a 
cold wind initially dominated by Poynting flux, the total 
energy flux per rest mass energy flux p can be expressed 
as 

„ ilrB.„Bd 



The second term of ( 12 ) c orresponds t o the Michel 
magnetization parameter cm ( |Michel 1969 )F] For a criti- 
cal solution in a monopole field geometry where the fast 
magnetosonic velocity is reached at infinity (x-p — > 00), 
the terminal Lorentz-factor becomes 



T(x F ) = /i 



1/3 



(14) 



r + E- 



Poynting 



Di fferent derivations o f this fundamental result are given 
by |Camenziiid| ( Il986[ ); |Tchekhovskoy et al.| ( |2009[ ). For 
small pert urbations fro m the monopole fie ld geometry 
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AirpiipC 2 

Foll owing the asymptotic relations by |Camenzind| 
(1986) we have v$ —> 0(x >• 1), and hence S2r ~ 
—B ! pB~ 1 Vp can be used to eliminate the toroidal field 
from equation [TT] With r> 1 we can write 

_ r _ n 2 r 2 B p 
Ankc 3 ' 

where /1, k = pu p /B p , and fl are conserved qua ntities 
along the stationary streamline (for details see e.g. Porth 

2 



Begelman fc Li| fl!994| ); |Beskin et ah] fl!998| ) could show 
that xf can be crossed at a finite distance, where again 



( 14) is satisfied. This genera l scaling was also found by 
5 Wat & Camenzindl (|1996D; |Fenat & Greinerl (pOOlT; 



Where p denotes the co-moving density, u v 



poloidal part ot the lour-velocity and h signifies the specific en- 
thalpy defined through the equation of state. 



Fenat fc (Juyeaj ( |2UU4p for coflimating relativistic jets. 
Our jet solutions quickly accelerate to the fast mag- 
netosonic point ana, aespite the aeparture from the 
monopolar shape, follow Michel's scaling there remark- 
ably well. As illustratea in figure [J] (right panel) , the 
aeviation from the expectation of pP z is less than 5%. 

Insight into the ongoing acceleration process can be 
gainea by an analysis of the trans-fieia force equilibrium 

3 Sometimes denoted as "field line bunching" in the recent liter- 
ature. 

4 Where we added the subscript "M" in order to avoid confusion 
with the parameter <r = cr M /T defined previously. 
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Figure 2. As figure [TJ for model 2h. 
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Figure 3. As figure [TJ for model 3h 

as performe d for example by |Chiueh et al.| (1991); Via- 
hakis (2004). The asymptotic reiativistic force balance 
can conveniently be decomposed into "curvature" , "elec- 
tromagnetic" and "centrifugal" contributions. Depend- 
ing on the dominating terms, at lea st two regimes are 
possi ble (see also the discussion by Komissarov et al. 
2009 1 : When the curvature term is negligible, the equilib- 
num is maintained by balancing of the centrifugal force 
with the electromagnetic contribution. This constitutes 
the first or linear acceleration regime. The transition to 
the second regime occurs when field line tension begins 
to dominate over the centrifugal force, maintaining the 
equilibrium between purely electro-magnetic forces. The 
occurrence of curvature in the force equilibrium leads to 
a tight correlation between collimation and acceleration 
since the tension force also becomes the governing accel- 
erating forc^] 

5 The latter was demonstrated using the parallel field force- 
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As far as a stationary state is reached, we find that the 
acceleration is well described by the linear acceleration 
regime far, or 



(15) 



et al 



as suggested for the initi al acceleration of rota t ing flows 

^ ( |1998|jN^aya^i 
and |Komissarov 



by various authors (e.g. Beskin et al. 



et al 



2007), Tchckhovskoy et al. (2008 



). This corresponds to 



B% » B% - E 2 



(16) 

which is satisfied for the most part of the flow in our 
simulation domain. Figure [5] (left panel) shows B 2 and 
B 2 , — E 2 for a sample field line in the fast jet. We 



balance i n appl ication to reiativistic disk wind simulations by |Porth| 
|fc Fendt| | |2010| l. 
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Figure 4. Acceleration along field lines rf p S {2, 4, 6, 8} shown re- 
spectively as {solid, dotted, dashed, dash-dotted} lines in model 2h. 



Left: The 



B p r function of the expanding flux tube normal- 



ized by footpoint value against radius in terms of light cylinder radii 
r/r\ c . Alfven (•) and fast (A) critical point transitions are marked 
accordingly. Right: Total energy flux ratio fi (top) and Lorentz 
factor at the fast point Fp compared to the expected value of fj. 1 ^ 3 
for various field line footpoints (bottom). We find Michel's scaling 
to be satisfied within 5%. 




1.00 10.00 100.00 

Figure 5. Characterization of the acceleration in model 2h. Left: 
Comparison of the field strengths for rf p = 2 along the flux tube. 
We find Bp > — E 2 for the most part of the domain yielding 
the linear acceleration regime. Also shown are power-law fits to the 
super-fast regime (thin solid lines). Right: The quantity rtan% 
along the same field lines of figure [4] Efficient acceleration in the 
power-law regime would yield rtan% ~ const. Alfven (•) and 
fast (A) critical point transitions are marked accordingly. 



find that the critical field strengths are fairly well ap- 
proximated by power-laws in the asymptotic super fast- 
magnetosonic regime. For the particular case shown, we 



have B 2 



cx r 



-4.3 



and Bl - E 2 



oc r such that the 
flow experiences a transition to the second, or power- 
law acceleration regime where the inverse of relation |16| 
becomes true. 

For the power-law regime, a correlation between 
Lorentz factor T and half-opening angle of the jet Or, 



T tan fl ~ 1 



(17) 



was discovered by Komissarov et al. ( 2009|) in the context 
of ultra-relativistic gamma-ray bursts! Figure [5] (right 
panel) illustrates the run of T tan 9a along a set of field 
lines in our fiducial model. Compared to the suggestion 



of equation 17 our simulation setup shows efficient MHD 
self-collimation, but appears less efficient in terms of ac- 
celeration. 

We note that only when a substantial part of flow ac- 
celeration takes place in the power-law regime, relation 
[T7]will hold. Our AGN jet models are however collimated 
to ~ 1° and accelerated with efficiencies of 40% (r ~ 8) 
already in the linear regime. Even if the flow acceleration 



is followed indefinitely, T6\f ~ 1 can not be recovered as 
this would require terminal Lorentz factors of F^ > 60 
and thus violate energy conservation. 

It could be argued that the low acceleration efficiency 
is due to the loss of causal connection for the relativis- 
tic flow. In this case, the bunching of field-line can not 
be communicated across the jet anymore, thus stalling 
the acceleration process. This should in fact occur when 
the fast Mach-cone half opening-angle Omf — t/2\/ M/T 3 
doe s not comprise the jet axis, hence Omf < #fl (see 
also |Zakamska et al.|2008| |Komissarov et al.||2009| . We 
have checked this conj ect ure by comparing both angles 
and found our still moderate Lorentz factor, highly colli- 
mated, jet models to be in causal connection throughout 
the whole acceleration domain. 



3.2. Thermal spine acceleration 

In this work, the very inner jet spine is modeled as 
a thermal wind. An alternative approach would be to 
prescribe a Poynting domin ated flow originating in the 
Blandford & Znajek (1977) process. In this case, the 
toroidal field would be generated by the frame dragging 
in the black hole ergosphere below our computational do- 
main similar to the induction in the disk. However, our 
attempts to increase the central magnetization a by fur- 
ther decreasing the coronal density failed at the inability 
of the numerical scheme to handle the steep density gra- 
dients emerging at the boundary. Due to the vanishing 
toroid al field at the axis, also the |Blandford fc Znajek| 
( |1977[ ) mechanism is not able to pro vide acceleration of 
the axial region (see also McKinney||2006[ ) . 

In principle, it would be possible to convert the thermal 
enthalpy first into Poynting flux when the jet is expand- 
ing, and then back into kinetic energy via t he Lorentz 



force as observed by Komissarov et al. (20091. However 



as we see in figures [T] to [5] (right top panels), this does 
in fact not occur in our simulations since the Poynting 
flux is approximately conserved along the inner flux lines 
that show little expansion. It is the magnetic field dis- 
tribution and the collimated structure of the outer (disk- 
jet) component which merely provide the shape of the 
trans-sonic nozzle for the thermal wind. We can thus 
understand the acceleration in the jet spine by using the 
relativistic Bernoulli equation, which we cast in the form 



h 2 [l + 2<p+(u/c) 2 ] = c 4 r 2 



const. 



(18) 



An order of magnitude estimate sufficiently far from the 
compact object yields hT ~ const. Using mass conser- 
vation Tpr 2 = const (v — > c) and a polytropic equation 
of state with the enthalpy h = c 2 + 7/(7 — l)p/p, we 
obtain a scaling relation F cx r^ 2+2 ^ 2 ^^ (p/p ^> c 2 ). 
For a relativistic polytropic index of 7 = 4/3 this results 
in r cx r. 

Applying a non-relativistic index of 7 = 5/3, the latter 
relation would yield T cx r 4 , however, the non-relativistic 
limit also implies p/p <C c 2 , and thus h — > c 2 (7 — > 5/3) 
and the acceleration ceases. 

Our simulations are perform ed employing the causal 
equation of st ate (obeying th e Taub ( 1948[ ) inequality) 
introduced by Mignone et al. (2005). Thus we obtain a 
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variable effective polytropic index 

dlnp (h — l)p/p 



7off 



d\np (h—l)p/p—l ^ 

between 4/3 and 5/3. Figure [6] (top panel) shows the 
effective polytropic index along a stream line / flux sur- 
face. In the sample stream line we find 7 er r to vary be- 




Figure 6. Thermal acceleration along the field line rf„ = 0.2 
in model 2h. Top: Effective polytropi c in dex 7 e ff along the flow. 
Bottom: Individual terms of equation |18| showing thermal energy 
conversion and the conservation of Too ~ 3.7. In this plot, we 
normalized to c=l. Alfven (•) and fast (A) critical point transitions 
are marked accordingly. 



tween 1.45 < 7 e jj < 1.65 as the plasma adiabatically 
cools from relativistic to non-relativistic temperatures. 
Thermal acceleration saturates for j e fi — > 5/3 as the en- 
thalpy approaches the specific rest mass energy c 2 (see 
also Fig. pi bottom panel). The maximum attainable 
Lorentz factor is given by the footpoint values at the 
sonic point to Too = /i (r + 20 o ) 1 ^ 2 and depends on the 
detailed modeling of the inner corona. In our approach 
the jet spine Lorentz factor is thus limited to values of 

roo <4. 

4. SYNCHROTRON RADIATION AND FARADAY 
ROTATION 

The numerical MHD simulations discussed above pro- 
vide an intrinsic dynamical model for the parsec-scale 
AGN core. In the following we will use this informa- 
tion - kinematics, magnetic field distribution, plasma 
density and temperature - to calculate consistent syn- 
chrotron emission maps. What is still missing for a fully 
self-consistent approach is the acceleration model for the 
highly relativistic particles which actually produce the 
synchrotron radiation. However, we have compared a 
few acceleration models and discuss differences in the 
ideal resolution synchrotron maps (see below). 

Radio observations of nearby AGN-cores show opti- 
cally thick and thin emissi on features with a high degree 
of Faraday rotation (e.g. Zavala & Taylor 2002). The 
nature of the Faraday rotation could either be internal, 
thus directly produced in the emitting volume, or due 



to an external Faraday sheet, possibly comprised of a 
magnetized disk- wind as ventured e.g. by |Broderick fc 
McKinney (2010), or an ambient jet cocoom On these 



scales, even with global VLBI experiments, the radio 
emission is barely resolved for most of the known sources. 
In order to confront the existing observations, we per- 
form linearly polarized synchrotron radiation transport 
in the relativistically moving gas, taking into account 
self-absorption and internal Faraday rotation. We apply 
beam averaging to examine the resolution dependence of 
the results. An illustration of our ray-tracing procedure 
with a rendering of an exemplary MHD solution of a col- 
limating jet is shown in Fig. [7] For a grid of lines of 



GAMMA 
7.837894 




Figure 7. Illustration of the ray-casting geometry on an exem- 
plary solution. Color-coding in slices (x = and z = —1000) and 
field-lines represents the bulk Lorentz-factor. An inclination angle 
i = corresponds to looking directly into the jet. 

sight, each corresponding to one pixel in the final im- 
age, we solve for the parameters of linear pola rization 
I = {P,I r ,U lr } as defined e.g. by |Pacholczyk| ( |l970[ ). 
This treatment provides the equivalent information as 
the Stokes parameters {I,Q,U}. Within the aforemen- 
tioned notation, the transport equation is a linear system 
of equations 



dl 



e a 1 



(20) 



where £ denotes the emissivity vector and A the opacity 
matrix, taking into account relativistic beaming, boost- 
ing and swing of the polarization as defined in Appendix 
[A") Faraday rotation of the re lativistically moving plasma 
Tias first been considered by Broderick & Loeb (2009b) 
and is directly incorporated into the previous relation via 
the observer system Faraday rotation angle 



dXF 
dl 



f(lt)n e D' 



(n-/3)-B' (21) 
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in cgs units, where e, m e ,n e denote the electron charge, 
mass and number density, z/, n the observed photon fre- 
quency and direction, D = (T (1 — n • /3)) 1 the Doppler 
factor and B' = B/T + T (/3B) • j3 the co- moving field. 
The dimensionless function f('jt) takes into account that 
for high plasma tempe ratures the n atural wave modes 
do not remain circular ( Melrose|1997 ) , suppressing Fara- 
day rotation in favor of conver sion between linear and 
circular polarization . We follow Huang et al. (2008) and 
Shcherbakov ( 2008 ) in defining 



/(7t) = 7t 



7t 



1 - 



ln7 t \ hi7t 



27t 



27t 



7t = 1 + 



ksTe 

m e c 2 
(22) 



in terms of the thermal electron Lorentz-factor j t to in- 
terpolate between the cold and relativistic limits. Es- 
pecially for the hot axial flow, Faraday rotation is thus 
substantially suppressed. Assuming an electron-proton 
plasma, the electron number density follows from the 
mass density of the simulations. We assume further that 
a small subset of these "thermal" electrons is acceler- 
ated to a power-law distribution and thus responsible for 
the non-thermal emission of synchrotron radiation. The 
modeli ng o f particle acceleration is detailed further in 
section 14.11 

A fraction of the Faraday rotation thus takes place al- 
ready in the emitting region of the relativistic jet, such 
that the radiation undergoes depolarization due to inter- 
nal Faraday rotation. In this case, the angular difference 
AXobs between the observed polarization angle Xobs and 
the (A — > 0) case can depart from the integral 



A* 



l B ■ dl, 



(23) 



which is customarily used in the diagnostics of jet obser- 
vations. For example, in a uniform optically thin medium 
with internal Faraday rotation, the value of A^obs is just 
half of relation |23| Non-uniform optically thin media will 
break the A 2 -law an d exhibit de polarization once A^obs 
exceeds 



Burn 1966). For optically thin me- 
AMaw rotation measures can be 



45° (e 
dia with Axobs < 45 
recovered also in the non-uniform case, however the ob- 
served rotation angle is always less than A^F 

We show the effect of internal Faraday rotation along 
an individual ray compared to the case with no Faraday 
rotation in Fig|8] In the emitting volume the polarization 
degree oscillates as expected for an optically thin medium 
with Faraday rotation. As a consequence, the observed 
polarization degree is lowered. Following the density and 
magnetic field strength, the differential Faraday depth 
dxe/dl decreases fast enough towards the observer, so 
we can be confident not to miss a substantial part of the 
Faraday screen in the ray-casting domain. 

To speed up the computation in cases of high optical 
or Faraday depths, we (i) limit the integration to t < 
100, and (ii) solve the polarized transport only for the 
last 200 internal Faraday rotations rp < 2007T. Both 
optimizations do not at all affect the resulting emission 
maps, as the observed radiation typically originates in 
the photosphere of r = 1, and only a few internal Faraday 
rotations suffice to depolarize the radiation in the models 
under consideration. 
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Figure 8. Raytracing for an individual lino of sight. Thick blue 
lines including Faraday rotation compared to a case where the lat- 
ter was neglected - illustrated by thin red lines. The upper panels 
show intensity (I tot ) and the Stokes parameters /' (dotted), I r 
(dashed) and U (second panel). Faraday depth dxp/dl [rad/ro], 
polarization degree ii and the optical depth r is shown in the sub- 
sequent panels. Internal Faraday rotation and the accompanying 
depolarization is observed in the emitting region near the x = 
plane (1 = 0). 

Once I is recovered, we obtain beam-averaged quanti- 
ties via the convolution 



<I)(x) 



J d 2 x'£(x-x')I(x') 



(24) 



with a Gaussian beam Q. The beam- averaged Stokes 
parameters are then used for mock observations provid- 
ing spectral indices, polarization maps, rotation measure 
maps, and spectra to be compared to the model param- 
eters. 

4.1. Particle acceleration recipes 

Within the MHD description of the jet plasma, knowl- 
edge about the relativistic particle distribution, which 
is needed as input for the synchrotron emission model, 
is not available. To recover the information from the 
velocity-space averaged quantities of MHD, we have to 
rely on f urther assump t ions. To mention other ap- 



proaches, Mimica et al. (20091 were able to follow the 
spectral evolution of an ensemble of relativistic particles 
embedded in a hydrodynamic jet simulation. Their treat- 
ment includes synchrotron losses, assuming a power-law 
seed distribution derived from the gas thermal pressure 
and density at the jet inlet. 

Alas, for our purposes a consistent prescription for in- 
situ acceleration and cooling would be required - which 
seems unfeasible at the time. We therefore take a 
step back and assume that relativistic electrons are dis- 
tributed following a global power law with index p as 
dn c = N Q E- p dE for E x < E < E u where N signifies the 
overall normalization of the distribution and E\, E u de- 
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note the lower and upper cutoffs. The optically thin flux 
density for the synchrotron process then reads S v oc v~ a 
with a = (p — l)/2. Optically thick regions radiate ac- 
cording to the source function S u = e„/ k v oc v 25 . 

This choice of particle distribution is justified by ob- 
servations as well as theoretical expectations for the par- 
ticle acceleration. The major physical mechanisms ca- 
pable of producing non-thermal relativistic electrons are 
(inte rnal) shock accel eration of relativistic seed electrons 
(e.g. Kirk et aL|2000[) and MH D processes like magnetic 
recon nection (TL yubarsky 2005]) or hydromagnetic turbu- 
lence (Kulsrud & Ferrari 19 7T[). Considering diff e rentia l 
rotation in relativistic jets |Kieger fc Mannhe im (2002); 
Rieger fe Duffy] ( |2004| ) ; \k\oy fc Mimical ( |2008| ) suggested 
particle acceleration by shear or centrifugal effects . 

In addition to Ei and E u , the normalization Nq de- 
pends highly on the mechanism under consideration. A 
straight-forward recipe is to connect the particle energy 
to the overall mass density (similar to Gracia et al.|2009 1 



N E- p dE, 



(25) 



where an ionic plasma consisting of equal amounts of 
protons and relativistic electrons is assumed. Thus, all 
available electrons are distributed following this relation 
and the Faraday effect is maintained by the "equivalent 



densi ty of cold electrons" oc E~ 2 (e.g. |Jones fe Odell 
19771 in contrast to relation (21 ) where we assumed that 
the most part of ele ctro ns is non-relativistic. 

An alternative to ( 25 1 is to specify the integral particle 
energy density and thus the first moment of the distri- 
bution funct ion. For their leptonic jet models, |Zakamska| 



et al. (2008) have assumed that the internal energy is 



carried by relativistic particles, hence the relation 



3p 



N E 1 ~ p dE 



(26) 



can be used to provide Nq from the gas pressure resulting 
from the simulations. 

In contrast to relativistic shock acceleration where the 
energy reservoir for the particles is the bulk kinetic en- 
ergy of the flow, MHD processes directly tap into the 
co-moving magnetic energy density, and can effectively 
accelerate the particles up to equipartition. Accordingly, 
for the equipartition fraction £b, 
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n e (E) 



N E 1 ~ p dE 



(27) 



E, 



is customarily used to estimate jet m agnetic field 
Konigl 1979) or vice ve rsa (ILyutikov et al.|2" 



fc 



strength from the observe d emission (e.g. IBlandford 
Konigl 1979) or vice versa ( |Lyutikov et al.[ 
ick fc McKinney||2010 ). To obtain peak fluxes m the Jy 



2005 



Broder- 



range, we have applied 6b = 0.1 for our fiducial model. 
For the following discussion we have adopted a = 0.5 
(p = 2), E u = 10 6 i?i and specified E\ = j t m e c 2 for appli- 



cation with re lation (25). For a = 0.5 and applying the 
recipes p6p7| , only the cutoff energy ratio E u /E\ enters 
logarithmically into the determination of Nq. Within 
these assumptions, the influence of the cutoff values on 
resulting jet radiation is marginal. The magnitude of E\ 
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Figure 9. Ideal resolution logarithmic I v / 'J^max maps for model 
2h at u = 43 GHz, i = 30° using various tracers for the relativistic 
particles: Density (left), thermal pressure (center) and magnetic 
energy density (right). Linear (e) polarization vectors are overlaid 
as white sticks. The (x-) scale is given in terms of ray-casting 
footpoint and corresponds to a physical extent of 1200rs. 

then merely determines the number density of relativis- 
tic particles, to be chosen consistent with the number of 
particles available for acceleration. In this first study, 
we neglect the spectral changes introduced by the cutoff 
energies as this would require a more detailed modeling 
of the particle content which is beyond the scope of the 
current paper. 

The observed morphology of the intensity maps is 
mainly given by the various prescriptions of iVo men- 
tioned above. In the following we briefly compare the 
resulting radio maps. 

4.2. Radio maps for different particle acceleration 
models 

Figure [9] shows ideal resolution maps for the aforemen- 
tioned particle acceleration tracers. For the sake of com- 
parison, Faraday rotation is neglected and with 43GHz 
we choose a high radio frequency to penetrate through 
the opaque jet base. 

All tracers show an almost identical polarization struc- 
ture, and highlight a thin "needle" owing to the cylindri- 
cally collimated axial flow with high density, and high 
magnetic and thermal pressure. The axial flow is slower 
than the Poynting dominated disk wind which is de- 
beamed and, thus, not visible at this inclination. Since 



the axial flow features /3||B (cf. A17), the resulting e po- 
larization vector reduces to the classical case, and points 
in direction perpendicular to the projected vertical field 
of the axial spine. In the case of the density tracer, the 
emission becomes optically thick, as indicated by the 
r = 1 contour. Correspondingly, the polarization de- 
gree is lowered and the direction of the spine polarization 
turns inside the r = 1 surface. Depending on the radial 
density and pressure profiles p oc r~ 3 / 2 , B 2 oc r~ 2 and 
the pressure distribution in the disk corona p oc r~ 5 / 2 , 
the emission at the base of the jet is more or less extended 
and dominates the flux in all three cases. Relativistic 
beaming cannot overcome the energy density which is 
present in the disc corona, and therefore necessitates a 
more elaborate modeling of the accelerated particles in 
the jet. This will be provided in section [5] 
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4.3. Relativistic swing and beaming 
In optically thin, non-relativistic synchrotron sources, 
the observed b polarization vector directly corresponds 
to the projected magnetic field direction of the emitting 
region and thus carries geometric information about the 
jet. This allows us to interpret parallel e vectors in terms 
of toroidal fields, w hile perpendi cular e vectors indicate 
a poloidal field (e.g. Laing 1980). Similar to all realistic 
models of MHD jet formation, our simulations feature a 
helical field structure that is tightly wound within the 
fast jet, but increasingly poloidal further out. Hence, 
the resulting polarization structure is that of the telltale 
spine and sheath geometry - across the jet, the polariza- 
tion e direction flips from being perpendicular to parallel 
and eventually returns to a perpendicular orientation. 

Due to aberrat ion and the accompanying swing of the 
polarization (e.g. Blandford fc Konigl|1979| , an interpre- 
tation in terms of pure geometrical effects is not longer 
applicable in flows with relativistic velocities, instead 
a kinematic jet model is required. For cylindric al (i.e. 



(z, 0)-symmetric) relativistic jets, Pariev et al. (20031 
have demonstrated how the optically thin polarization 
follows a strictly bimodal distribution, since the inclined 
polarization vector from the front of each annulus can- 
cels with the corresponding polarization vector from the 
back side. This remains also valid for differentially ro- 
tating jets. 

For the case of a collimating and accelerating jet as 
shown here, we loosen the constraint of cylindrical sym- 
metry to mere axisymmetry in the (f> direction. Addition- 
ally, our simulations feature a non-constant pitch of the 
magnetic field and a small degree of rotation. Together, 
this results in inclined polarization vectors which deviate 
from t he s trict bimodality observed in (z, (/>)-symmetry. 
Figure 10 shows the optically thin polarization structure 
in the presence of relativistic effects (left panel), and 
in absence thereof (right panel). To produce the non- 
relativistic map, we had simply set v = before conduct- 
ing the radiation transport. At the base of the outflow 
where the velocities are only mildly relativistic, the po- 
larization vectors are found predominantly perpendicular 
to the collimating poloidal magnetic field in both cases. 
Further downstream, the bimodal spine-and-sheath po- 
larization structure prevails as the jet dynamics becomes 
increasingly cylindrical. The figure clearly demonstrates 
how the relativistic swing skews the spine towards the 
approaching side of the jet. Note that the jet rotation is 
also apparent in the beamed asymmetric intensity con- 
tours. At high inclination i = 60°, the main emission 
from the high-speed jet is de-beamed and only the low- 
velocity "needle" of the thermal spine along the axis can 
be recognized. It is worthwhile to note that both in- 
tensity and polarization of axisymmetric non-relativistic 
synchrotron sources exhibit a p oint symmetry about the 
origin as illustrated in Fig. [TU] 

4.4. Pitch-angle dependence 
Several Authors ( Marscher et al.|2002 Lister & Homan 



2005 1 found indications tor a bimodal distribution ot the 
electric vector position angles (EVPA) of quasars and 
BL-Lac objects either aligned or perpendicular to the jet 
direction. It was also supposed that BL-Lacs tend to 
aligned e vectors and overall higher degree of polariza- 




Figure 10. Relativistic effect on optically thin polarization (run 
2h) at i = 60°. The polarization degree is indicated by the back- 
ground coloring, black contours show total intensity levels spaced 
by factors of two. Left: Including relativistic aberration. At this 
high inclination, the jet is de-beamed and rotation is apparent in 
the asymmetry of the intensity contours. Right: In the absence 
of relativistic effects, the polarization pattern is point symmetric 
about the origin and the intensity clearly promotes the jet. 

tion. In a re c ent 8 6 GHz polarimetric survey however, 
Agudo et al. (2010) found no such correlation in their 
fiat-radio-spectrum AGN sample, rather are their data 
consistent with an isotropic (mis-)alignment. Alignment 
is customarily attributed to oblique shocks, while perpen- 
dicular EVPAs are then interpreted in terms of a shearing 
of the magnetic field with the surrounding. Another pos- 
sible explanation for the bimodality is due to large-scale 
helical fields as shown in the previous section. By varying 
the emitting region within the collimating jet volume, we 
investigate to what degree the polarization still conveys 
the geometric information of the emitting region. As 
before, Faraday rotation is neglected and we introduce 
the co-moving pitch angle tan = B'.jB' in analogy to 



Lyutikov et al. (20051. In a non-rotating, cylindrical jet, 
the pitch-angle transforms as 

(28) 



tan ^ = T tan 



and hence the co-moving fields appear much less t wist ed 
than their laboratory frame counterparts. Relation[28]to- 
gether with equation[l5]for the linear acceleration regime 
yields the co- moving fields 



1. 



(29) 



Substantially higher pitches are not realized within the 
simulations. An impression on the pitch angle distribu- 
tion throughout the jet can be obtained with the cuts 
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shown in the middle panels of figures [T] to [3] We restrict 
the emission to originate from B'^/Bp > {1,2}, corre- 
sponding to > {45°,63°},(# > {83j\86°}) where 
equipartition particle energy density (eq. 27 ) is assumed. 
This way, only the regions of the current driven jet con- 
tribute to the emission and no radiation is obse rved from 
the spurious axis where must vanish. Figure 11 shows 
the polarization for the two cases and various viewing an- 
gles. 

In the high pitch-angle case B'^/B'p > 2, the resulting 
EVPA become parallel for viewing angles i > 10° while 
for the standard case B'^/B' p > 1, this happens only at 
viewing angles i > 30°. Most structure is observed at 
moderate viewing angles where we can find the core po- 
larization perpendicular to the ridge line polarization for 
the c ase i = 30°, B'^/B^ > 1 as observed in some sources 
(e.g. |Pushkarev et al.||2005). Here we also see a spine 



and sheath polarization profile across the jet. In the 
adopted parametrization of the emission region, the spine 
and sheath structure is only observed for B'^jB' p ~ 1 as 
higher pitches tend to alignment and lower pitches tend 
to counter align. 

Beam depolarization is apparent in regions where the 
polarization turns and, consistent with most observa- 
tions, the degree of polarization increases towards the 
boundary of the jet. In the high-pitch left-right 
asymmetry is most significant with parallel vectors on 
the approaching side of the r otating jet and perpen dicu- 
lar ones at the receding one. |Lyutikov et al. ( 2005 ) have 
proposed that based on the asymmetric polarization sig- 
nal, the handedness of the magnetic field and thus the 
spin direction of the black hole / accretion disk can be 
inferred. Our results support this finding. Clearly, when 
increasing the pitch, the EVPAs turns from the perpen- 
dicular to the parallel direction as a general trend. How- 
ever, pitches B'^j B' p > 4 were not realized by our simula- 
tions, such that the emitting region would vanish. Thus 
for the models under consideration, when approaching 
the Blazar case, the intrinsic polarization always appears 
perpendicular. In the following we will mostly character- 
ize the emission region by B'^jB^ > 1 as this best selects 
the relativistic jet contribution. 



5. RADIATION IN THE JET MODELS 

We have now introduced all parameters to build mod- 
els for the synthetic observations. A radiation model 
comprises a MHD simulation run, a particular physical 
scaling given by M, , E and the specific parameters of the 
radiation transport "f^ee and i. To normalize the flux, 
a photometric distance D = 100 Mpc is assumed in all 
models. Table [2] summarizes the parameters adopted in 
this work. We will mostly report results for the fiducial 
model A and consult the other models only for compar- 
ison with the standard case. 

For the interpretation of the results, an understanding 
of the qualitative dependence of the observables on our 
parametrization is helpful and therefor discussed in the 
following. In terms of the physical scaling and equipar- 



tition fraction, we have 

e u oc e B E 7 ' A M. 7/2 i/~ Q (30) 
K V oc e B E 2 M~ 4 v - a -^l 2 (31) 



dx 



dl 



£ cx e 3 ' 2 m: 



(32) 



The observable quantities then become 

C in cx e v x / cx e B £ 7/4 M. 5/2 v~ a (33) 
Il hick cx e v /K u cx E 1 ^ Ml 12 v 2 - 5 (34) 
and the opacities follow to 



T v CX OL v X I cx £b E 2 M % 6 V 
Xf oc E 3/2 M, 2 . 



-a-2.5 



(35) 
(36) 

In reverse, based on the latter relations, spectrum and 
Faraday rotation measurements will allow us to constrain 
the physical parameters]^] 

We deliberately chose model B to feature ~ 32 times 
higher Faraday depth compared to model A by increasing 
the jet energy E, while maintaining a similar spectrum 
with the choice of £b (although the total flux is thus 
decreased by a factor of ~ 1.8). 

5.1. Core shift 

Lobanov ( 1998 ) showed how opacity effects in optically 
thick jet cores provide valuable information that can help 
to constrain the dynamical jet quantities. As the bulk of 
the radiation originates in the photosphere r = 1, the 
projected distance of the r = 1 surface will result in 
a specific core offset. Due to the frequency dependent 
opacity of a synchrotron self absorbing radio source, the 
measured core position varies systematically with the ob- 
serving frequency. In simple model jet where the mag- 
netic field and relativistic particle densi ty is modeled as 
B cx r~ m and N cx r~ n (Konigl 19811, this projected 
distance becomes 



r coro cx v- 1 '^ (37) 

where k r is a combination of the parameters m, n and 
the spectral index a 



(3 + 2a)m + 2n-2 
5 + 2a ' 



(38) 



In a conical jet, the (predominantly toroidal) magnetic 
field strength will follow m = 1 due to flux conservation. 
Conservation of relativistic particles sets for the number 
density n — 2. It is noteworthy that in this most reason- 
able case, the exponent k r becomes independent of the 
spectral index k T = 1. Although our jet models are col- 
limating and thus not exactly conical, qualitatively we 
expect fc r clo se to unity when equipartition particle en- 
ergy (eq. 27 ) and thus n = 2m is assumed. The observed 



core shift is illustrated in figure 12 for the fiducial model 
A. Fits of the k r exponent for all runs are shown in the 
right panel. 

6 In practice this is further complicated by the dependences in- 
troduced by the doppler factor which we omitted here as we will 
compare only dynamically identical models at a given inclination. 
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Figure 11. Polarizations for i £ {30°, 20°, 10°} (from left to right) emitted from regions with co-moving pitches B'^/B'p > 1 (above) and 
B'^/B'p > 2 (below). The polarization degree II43qh z is color-coded and /43GHz contours are shown. Contours are spaced by a factor of 2 
out to ~ 5 ■ 10 — 4 7 !/iPea k where the image is cropped. Spatial scale is given in milli arcseconds and a restoring beam with FWHM=0.05 mas 
was used. The right-hand panel shows polarization angles along the cuts along core (black) and jet (gray). 



Table 2 

Jet radiation models 



model ID 


run ID 


M.[Mq\ 


E[erg/s] 


CB 


tan*' 


i[deg] 
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-^Edd 
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2h 


10 9 


10 43 


0.1 


1 


20 


1.3 x 10" B 


13 


2h 


10 9 


10 44 


0.001 


0.5 


20 


1.3 x 10~ 5 


C 


Hi 


10 9 


10 43 


0.1 


1 


20 


5.5 x 10~ 7 


D 


3h 


10 9 


10 43 


0.1 


1 


20 


1.1 x 10~ 6 



Note. — Columns denote: Model ID; Simulation ID; Black hole 
emission; Standard inclination; Resulting jet mass loss rate in terms of 
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Figure 12. Left: Half maximum intensity contours and peak po- 
sitions illustrating the core shift in mas in model A for different 
frequencies: 344GHz (blue) to 15 GHz (maroon). The images are 
aligned with respect to the imaginary black hole line of sight lo- 
cated at the origin. Right: Fit of the core distances relative to the 
344 GHz peak for models C,AjD (top to bottom). The fit function 
reads Ar = A(v~ 1 / kr — 344 _1//fer ) with free parameters A and k T . 



mass; Total energy flux rate; Minimal co-moving pitch angle of 
the Eddington accretion rate (with a radiative efficiency of 0.1). 

The best fits are slightly steeper than the conical ex- 
pectation of k r = 1. Within the fitting error, models C 
and D are however still consistent with k T = 1 while the 
deviation in model A becomes significant. It is tempting 
to interpret this behavior as a consequence of jet colli- 
mation which increases the shift for larger distances from 
the origin. As the far regions are probed by lower fre- 
quencies, the distance law of a collimating jet is expected 
to systematically steepen towards the low frequency side. 
The quality of the fits for our model jets suggests that 
core shifts can provide a robust diagnostic of the AGN 
jet acceleration region. Unfortunately, core shift largely 
complicates the interpretation of th e Faraday rotation 
maps, as we will show in section [5H} 

5.2. Depolarization 

To understand the polarization signal of the simulated 
jets, let us briefly review some considerations on the po- 
larization degree in general (e.g. Pacholczyk 1970 ), in the 
presence of Faraday rotation (e.g. Burn 1966) and when 
an observing beam is used. 

In the simplified case of a uniform magnetic field and 
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for power-law electron distributions, the polarization de- 
gree reads 



n h = 



(p + i) 

(p + 7/3) 



0.69 (p = 2) 



for the optically thin regime and 
3 



n, 



(6p + 13) 



0.12 (p = 2) 



(39) 



(40) 



for optically thick radiation. Such high polarization de- 
grees are however rarely observed in AGN cores which 
go down to the percentage level. Clearly a mechanism 
for depola rization is req uired. 

Already Burn (1966) considered the admixture of a 
isotropic random held component B T to an otherwise 
ordered field Bo in the emitting region. He found the 
simple relation for depolarization relating the energies of 
the fields II r ~ n Q B^(B^ + B^)~ l where II r is now the 
reduced polarization due to the additional random com- 
ponent. In result, to obtain significant depolarization 
the random component needs to be comparable to the 
ordered field. The accompanying dissipation of energy 
would notably decrease the efficiency of the jet accelera- 
tion process, increasing the scales of jet acceleration pos- 
sibly beyond the parsec scale. The occurrence of turbu- 
lence and thus randomly oriented fields in the AGN core 
can serve as an explanation for various multifrequency 
observations as ventured by Marscher (2011). Under 



additional compression, even a completely randomized 
field structure can account for the high an d low observed 
polarization degrees as demonstrated by Laing (1980). 
However, the complex physics of turbulence within the 
jet can not be incorporated into MHD simulations of rel- 
ativistic jet formation at this time. H Our current MHD 
simulation models thus provide highly ordered near force- 
free fields around which relativistic electrons following an 
isotropic distribution are assumed to gyrate. We propose 
that the most promising site for finding such ordered (he- 
lical) fields is in fact the Poynting dominated regime of 
the jet acceleration region that is simulated here. 

Under the influence of Faraday rotation within the 
emitting volume, the polarization degree will also depend 
on the Faraday opacity /3p = d\F/ds. 
Considering first optically thick radiation where the po- 
larization degree is governed by the ratio of ordinary to 
Faraday opacity 6 = k u /0f. Specifically, with k„ cx 



,-p/2-2 and p F 



it becomes 5 * 



,-p/2 



(P 



2). Once the photon mean free path is smaller than 
the correlation length of the field and the mean rota- 
tion length /Jp 1 in the low frequency limit (5 oo), we 
expect n to approach the uniform magnetic field case 
ETj. Accordingly, for a small value of 6, the radiation is 
depolarized II -> (5 -> 0). 

For optically thin radiation, the impact of Faraday ro- 
tation can be parametrized by rj = fi-ps describing the 
angular change during the emission length s. The polar- 
ization degree will then decrease and oscillate according 
to n = EEh| sin 77/77 1 . This is known as differential Faraday 
rotation and its influence on n is shown in figure [8] along 

7 Note on the other hand that simulations featuring tur bulent 
slow disk winds ser ving a s Faraday screen were presented by |Brod^] 
|erick fc McKinney | ( [2010] > . 



a line of sight in the simulation. Since 77 oc v~ 2 , internal 
Faraday rotation has vanishing influence also in the high 
frequency limit and n approaches the maximal polariza- 
tion degree In a non-uniform field, changes in the 
emitting geometry lower the maximal polarization degree 
and Uh can only be assumed at the edges of the emission 
region. The direction of preferred emission perpendicu- 
lar to the projected magnetic field is also the direction 
of dominant absorption, such that in the uniform field 
case the optically thick polarization direction is flipped 
by 90° with respect to the optically thin polarization. 

Additional depolarizing effects occur when an extended 
source is observed with finite resolution. The radiation 
is depolarized when the beam encompasses (1) intrin- 
sic changes in the emitting geometry, (2) optical depth 
transitions leading to 90° flips and (3) varying Faraday 
depths, know n as (internal or external) Faraday disper- 
sion (see also Gardner & Whiteoak 1966; Sokoloff et al. 



1998) 



For completeness, we should also mention depolariza- 
tion via "blending" - contamination with an unpolarized 
(thermal) component - for example radiation from the 
torus in the infrared. This occurs when the intensity of 
the contaminant becomes a notable fraction of the total 
intensity and thus imprints on the spectral energy dis- 
tribution as well. Naturally, our results are only valid as 
long as the jet-synchrotron radiation is dominating the 
total flux. 

5.2.1. Low Faraday rotation case 

Spectrum and core polarization degree for model A 
is shown in figure |13| The observable core polarization 
degree is defined as 



J dn (i v (Q)) n((i„(n)» 

/ dVLl v (VL) 



(41) 



where the quantities under the integral are themselves 
subject to beam convolution. In absence of the lat- 
ter, the averaged polarization degree increases monoton- 
ically from the expected value of ~ 0.1 in the low fre- 
quency range to 0.4 in the optically thin case. Due to 
the varying emitting geometry, the theoretical maximum 
of Hh = 0.69 is not realized. The influence of differen- 
tial Faraday rotation seen in the ideal resolution case is in 
fact small. Depolarization occurs through ordinary beam 
depolarization and through Faraday dispersion, once the 
angle between the (unresolved) emitted polarization di- 
rection and the Faraday rotated polarization bec ome s 
larger than ~ 45° (compare with lower panel of Fig. 13 ). 
We note that also for the unrotated polarization direc- 
tion, a flip of ~ 90° between the thick and thin case 
is not observed. This is also expected, since the pho- 
tosphere probes various pitch angles as the frequency is 
decreased and so the simple uniform field case does not 
apply. 

With increasing beam size, the polarization degree ul- 
timately converges to the unresolved case. The conver- 
gence is faster for the optically thin regime where the in- 
trinsic emission is less extended and thus quickly masked 
by the beam. 

5.2.2. High Faraday rotation case 
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Figure 13. Spectrum, core polarization degree and core polar- 
ization direction in model A showing depolarization due to beam- 
and Faraday effects. The curves in the middle panel are shown 
for increasing beam-FWHM in mas as indicated and converge to- 
wards the unresolved case. The unresolved polarization direction 
is shown in the lower panel. Thick blue lines take Faraday rotation 
into account and thin red lines are calculated in its absence. 

To observe the effect of depolarization due to internal 
Faraday rotation, we perform the same analysis as be- 
fore for the high Faraday rotation model B. Figure 14 
shows the resulting quantities. As anticipated, a simi- 
lar spectrum is obtained but the polarization degree and 
direction behave differently. Even when observed with 
ideal resolution, the polarization degrees of the two cases 
(Faraday active vs. neglected Faraday rotation) sepa- 
rate clearly as a result of internal depolarization. Also 
the unresolved polarization angles separate at higher fre- 
quencies. Due to multiple rotations, the polarization an- 
gle appears to fluctuate below observing frequencies of 
100 GHz. In principle, the ideal resolution polarization 
degree is expected to rise again for lower frequencies as 
S — > co. However, this did not yet occur at frequencies 
above 4 GHz that were under investigation. 

5.3. Rotation measure 

A helical field geometry is generally perceived to pro- 
mote transversal rotation measure (RM) gradients owing 
to its toroidal field co mponent . First evidence for RM 



gradients was fo und by Asada et al. (2002) and | Zavala 
fc Taylor| ( |2005[ ) in the jet of 3C273. In several unre- 
solved radio cores, A 2 law RMs have been detected and 
are found to fol l ow a m onotonic profile for example by 
Gabuzda et aTl pOOl) ); |Q'Sullivan fc Gabuzdaj ( |2009| ; 
Croke et alJlg m OTT 

Taylor fe Zavala| (|2010 1 point out the observational re- 
quirements of a RM gradient detection as follows: 1. 




Figure 14. 
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As figure |13| but for model B with high Faraday 



At least three resolution elements across the jet. 2. A 
change in the RM by at least three times the typical er- 
ror. 3. An optically thin synchrotron spectrum at the 
location of the gradient. 4- A monotonically smooth 
change in the RM from side to side (within the errors). 

In the following paragraphs we will touch up on each 
of the aforementioned points. Due the flip between the 
optically thick and thin polarization direction, the mea- 
surement of A 2 -law RM around the spectral peak require 
extra caution. Using sufficiently small spacings in AA 2 
to recover mr rotations, we can fit the rotation measure 
law 



X(A 2 ) = Xo + RMA 2 



(42) 



to the optically thick and thin cases, where xo now de- 
notes the effective angle of emission. The fits are shown 
for a particular line of sight in figure |15| Here, xo dif- 
fers by almost 90°, whereas RM is ofcomparable size. 
However, the latter two findings are not necessarily true 
for all lines of sight, since the optically thin photons can 
originate in higher Faraday depths and different emitting 
geometry, as mentioned previously. We stress that with 
ideal resolution, consistent A 2 -laws are found both for 
optically thick and thin photons. Alas, if taken together, 
we would not be able to fit a linear function for the whole 
range. 

The two-dimensional RM maps shown in figure [16] 
demonstrate the effect of beam averaging on the low- 
and high frequency regime. In the high frequency case, 
the intensity is strongly peaked close to the central object 
where the Faraday depth is highest, leading to steep ra- 
dial gradients also in the rotation measure. As the emis- 
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Figure 15. Determination of the rotation measure for an ideal 
resolution line of sight. The value of RM (indicated in rad/m 2 ) 
necessitates multiple rotations by n for the high A 2 case. In the 
inlay, a different fit for the optically thin regime is shown. Between 
the two cases, the effective angle of emission \0 (indicated in deg) 
is rotated by 83°. 
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Figure 16. Rotation measure maps (i = 20°) in the optically 
thick (top) and in the optically thin regime (bottom). White sticks 
indicate the direction of effective emission, xo- In the top plot inlay, 
only the 22 — 86GHz region is fitted while the additional points 
illustrate the A 2 -law breaking due to the core shift (see text). 

sion is more extended, the RM is lower and smoother 
in the low frequency case. We find an interesting re- 
lation between the core shift and the rotation measure: 
As the photosphere moves outwards with decreasing fre- 
quency, the flux tends to originate further away from the 
central object leaving systematically less Faradav active 
material between the observer and the source. |3 After 

8 As a limitation of our direct ray tracing method, the photo- 
sphere can in principle shift to the outer boundary (the "lid") of 



beam convolution, the polarization angle is biased to- 
wards the outward moving photosphere, resulting in a 
shallower rotation measure for lower frequencies which 
ultimately breaks the A 2 -law that is valid for each ideal 
resolution line of sight, (e.g. figure 15). 

Once the core shift distance is large compared to the 
scale of typical changes in the Faraday depth or emission 
angle, a A 2 -law will not be observed when the observa- 
tions are aligned according to the core position. Even ab- 
solute positioning as done herqjcan successfully reestab- 
lish the A 2 -law only when the beam size is smaller than 
the core shift. In practice, unresolved core shifts could 
well be the origin of optically thick A 2 -law breakers. 

5.4. Resolution, mm- VLBI 

In order to observe the helical fields of the jet acceler- 
ation region via the associated rotation measure, beam 
sizes able to resolve the dynamics are essential. In ad- 
dition, the frequency must be sufficiently high to peer 
through the self-absorption barrier. With the advance of 
global mm- VLBI experiments, substantial progress will 
be made on both of these fronts. The theoretical reso- 
lution of a 10 4 km mm-observatory (at 300GHz) evalu- 
ates to lO^as, similar to the resolution of space- VLBI at 
86GHz. Corresponding to a physical scale of ~ 60 rg, our 
reference object is thus entering the regime of interest for 
RM studies. 

At present, the record holders in terms of physical and 
angular resolution are the 1.3mm observations of the ra- 
dio source Sgr A* ne ar the galactic ce nter black hole 



that w ere reported by Doeleman (2008) and Fish et al 
(2011 1. Coherent structures on scales less than 45/xas 
or ~ 4rs could already be detecteq 10 | and sub-mm rota- 
tion measure magnitudes in excess of4 x 10 5 rad/m 2 were 



discovered and confir med by Macquart et al. (|2006|) and 
Marrone et al] ( |2007| ). 



For t he prominent case of M87, [Broderick fc Loeb] 



(|2009ah elaborate that the inner disk and black hole sii- 
houette at 5 r$ could be observed with (sub) mm- VLBI. 
In M87, high frequency radio observa tions are already 



pushing towards the horizon scale (e.g. Krichbaum et al. 



2006 Ly et al. 



tion signal is s 



2007), only the imp ortant core pTj lariza- 



ill inconclusive (e.g. Walker et al.| 2008) 
Rotation measures at core distances ot ss 20mas vary be 
tween —5000 a nd 10 4 rad/m 2 dependin g on the location 



in the jet (e.g. Zavala fc Taylor] 2002). It is tempting 
to extrapolate these values to the /xas scale, assuming 
the observed RM values are non-local enhancements due 
to cold electron over de nsities. Following this argu ment 
that was put forward by Broderick & Loeb (2009a), the 



resulting core RM's would be on the order of 10 s rad/m 2 , 
enough to account for the low observed polarization de- 
gree via Faraday depolarization. 

In the mm wavelength range, large rotation measures 
are needed to produce detectable deviations in the po- 
larization angle due to the decreasing coverage of A 2 

our domain, such that the very optically thick regime below 8GHz 
can not reliably be probed. 

9 Our images are aligned with respect to the imaginary black 
hole line of sight. In practice, optically thin features should be 
used as indicators for an absolute alignment. 

10 Given several Jansky flux in the mm-range, imaging the black 
hole shadow in Sgr A* be comes a mere problem of visibility in the 
southern hemisphere (e.g. |Miyoshi et al.|2007| . 
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space. Typical calibration errors of ~ 1° require RM > 
18 x 10 3 rad/m 2 for a 3<r detection between 172 and 688 
GHz and RM > 4.3 x 10 3 rad/m 2 when 86 GHz obser- 
vations are added. The steep spectrum synchrotron flux 
of the jet rapidly declines when higher frequencies are 
considered, however, opacity and Faraday depth also de- 
crease to yield a higher c ontr ibution of polarized flux 
from the core (e.g. Figure 13). Also the deviations in- 
troduced by the core-shift as described in the previous 
section will pose lesser problems in (sub-) mm observa- 
tions. 

5.4.1. sub-mm Rotation measure maps 

To obtain detectable deflections of the polarization an- 
gle also at short wavelengths, we now focus on the h igh 
Faraday depth model B. As discussed in section 5.2.2[ at 



frequencies beyond 172 GHz, the unresolved polarization 
vector exhibits changes below 45° and internal depolar- 
ization is not observed (see also Figure 14), such that 



we find consistent A 2 -law rotation measures in the mm 
wavelength range. For an increasing beam FWHM from 
6.25/ias to lOO^as, sub-mm RM maps for the optically 
thin radiation are shown in figure |17| The correspond- 
ing physical resolution is 31rs — 500rs. Steep gradients 
of RM across the jet axis and "spine and sheath" po- 
larization structures are observable down to a resolution 
of 125 rg, below which most information is destroyed by 
the beam. With increasing beam size, only the central 
Faraday pit remains detectable. Here the core rotation 
measure reaches values as high as 10 6 radm~ 2 . 

We show the transversal cuts along the core and jet in 
figure 18 At a beam size of 500rs and observing frequen- 
cies between 172 — 688GHz the transversal RM gradients 
of the jet fall below the assumed detection limit of 3° 
and become consistent with a constant. Interestingly, 
beam convolution decreases the magnitude of RM not 
only in the cuts exhibiting a sign change, but also for 
the cuts across the intensity peak. We note that the in- 
trinsic jet RM profiles are non-monotonic. Only when 
under-resolved, the RM features the monotonic profiles 
"from side to side" that are typically observed. 



5.5. Viewing angle 

To investigate the viewing angle dependence on the 
observations we show radio observabl es f or different in- 
clinations from 0.01° to 40° in figure [l9j With a reso- 
lution of 50/ias in the 43-86 GHz frequency range these 
images preview the next generation space VLBI exper- 
iments. We observe asy mmetric features in the spec - 
tral index as proposed by Clausen-Brown et al. (2011). 
Due to the core shift, the maximum of the spectral in- 
dex appears "behind" the low frequency intensity peak. 
When seen right down the jet, the polarization vectors 
are radially symmetric and show an inclination about the 
radial direction. With increasing viewing angle, the pre- 
dominating polarization direction with respect to the jet 
flips from perpendicular to parallel orientation. Also the 
transversal polarization stru cture exhibits asymmetries 
as mentioned in section |4.4| To produce rotation mea- 
sure maps, we fitted the A z law to observations at 43, 85 
and 86 GHz. In order to avoid the problems due to op- 
tical depth effects (see section 5.3) we exclude regions 
of spectral index between < a < 2 in the maps. The 



strong feature in the 20° RM map at (0, 0.15)mas is most 
likely an artifact from the finite ray-casting domain and 
corresponds to the region where the axis pierces through 
the "lid" of the domain. For i — 40° this region is ex- 
cluded from the RM map and for i — 0.01 the problem 
does not arise. At high inclinations, we observe steep 
RM gradients that coincide with the spine-sheath flip of 
polarization. 

As such high resolution data are not readily avail- 
able, we also quantify the integral values that should 
reflect unresolved core properties. The viewing angle de- 
pendence of flux, spectral index, unresolved polarization 
degree, polariz ation direction and rotation measure are 
shown in figure 20 The flux peaks at viewing angles of 
roughly 5° and not when looking directly down the jet 
as would be the case for simplified cylindrical flows with 
toroidal fields. This reflects the fact that our jet solutions 
are not perfectly collimatcd and exhibit a small degree 
of rotation such that the Doppler factor attains its max- 
imum value of D ~ 15 at i ~ 5°. The dominating polar- 
ization direction flips from perpendicular (0°) to parallel 
(90°) for viewing angles i > 25°. At this flip, the po- 
larization degree shows a local minimum through beam 
depolarizati on. In the framework of the unified model 
proposed by Urry & Padovani (1995), this suggests that 
the core polarization direction in the radio galaxy case 
(viewed at high inclination) should be clearly distinct 
from the blazar case (looking down the jet). For i = 0°, 
the polarization degree approaches zero and its direction 
fluctuates between 0° and 90° due to the axial symmetry. 

5.6. Towards modeling actual observations 

Individual sources are modeled by constraining all free 
parameters (tables [I] and [2j with observations of spectra, 
core shift, polarization ana rotation measures. Taken to- 
gether, we obtain seven interrelated parameters in our 
radiation models. 

Two of these are related to jet dynamics, parametrizing 
the energy partitioning with /3\ and the current distri- 
bution with the parameter s. We find that the current 
distribution has little influence on the observables and 
hence this dependence can eventually be dropped, which 
reduces the number of free parameters but at the same 
time looses predictive power. 

Two parameters are related to the physical scaling of 
the simulations, namely the total energy flux E and the 
black hole mass M # . With the scale free nature of the 
underlying MHD, a single dynamical simulation can be 
used to construct a multitude of objects for arbitrary to- 
tal energy flux and black hole mass. Thus in principle a 
library of physical jet models can be constructed from a 
few dynamicall y di stinct simulations (see also: Physical 
scaling, section 2.41). 



Compared to the simulations, the raytracing is fast and 
can be used to vary the remaining parameters £b > and 
i. 

An alternative to modeling individual sources is to 
compare the statistical properties of a set o f models with 
a sam ple of AGN cores as observed e.g. by |Agudo et al.| 

poioF . 

Such an effort is beyond the scope of the current paper 
and we leave this open for future investigati on. However, 
already with the acquired data (e.g. figure 20), our sim- 
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Figure 17. RM maps in the optically thin wavelength range (344 — 68 
6.25/xas to 100/^as doubling the beam size for each image. The effective 




8 GHz) of model B when observed with decreasing resolution from 
emission angle \0 is indicated by white sticks. 

jet dynamical structure. 



0.10 



Figure 18. RM cuts along the paths indicated in figure |17| 
for various beam sizes. The physical resolution of the cuts is 
62.5rg — 500rg as indicated and 3° detection limits for 86 — 688 GHz 
and 172 — 688 GHz observations are shaded grey in the top panel. 
Curves are crop ped where the intensity falls below 5 X 10~ 4 7 l/ p ea i c 
as in figure [l7| 

ulations strongly suggest a bimodal distribution of the 
polarization angle. 

6. SUMMARY AND CONCLUSIONS 

We have performed axisymmetric, special relativistic 
MHD simulations of jet acceleration and collimation. 
The resulting dynamical variables are applied to calcu- 
late polarized synchrotron radiation transport in post- 
processing, providing emission maps consistent with the 



6.1. Jet acceleration 

Our jets are realistically modeled to consist of two 
components: an inner thermal spine assumed to origi- 
nate in the a black hole corona, and a surrounding self- 
collimating disk jet driven by Poynting flux. 

We follow the flow acceleration for more than 3000 
Schwarzschild radii reaching Lorentz factors in the disk 
jet of r ~ 8 within the AGN "blazar zone" where we 
calculate the synchrotron emission maps. In application 
to a 10 9 Af Q black hole this translates to a distance of 
0.3pc. 

Although the Poynting dominated jet flow becomes 
super fast-magnetosonic within the domain, it has not 
yet reached equipartition between Poynting and kinetic 
energy - jet acceleration is still ongoing. According to 
the available energy budget, in the case of high energy 
disk jets terminal Lorentz factors of Too ~ 20 would be 
aquired asymptotically. At the fast magnetosonic point 

1 /3 

we find the Michel scaling T(xf) = Tob to be satisfied 
within 5%. We find that the jet acceleration up to a 
distance z ~ 3000rs is well described by the linear re- 
lati on r oc r as proposed by Tchekhovskoy et al. (2008) 



and Komissarov et al.| ( |2009[ ). We do however not repro 
duce the tight coupling of acceleration and collimation 
r tan 8r ~ 1 observed in the latter communications but 
instead find T tan 9^ to monotonically decrease along the 
flow. The fast jet component in all models considered 
collimates to half-opening angles of ~ 1°. We find the 
causal connection within the flow - its ability to commu- 
nicate with the axis via fast magneto-sonic waves - to be 
well maintained in our simulations. 

Also the thermal spine acceleration is shown to be ef- 
ficient with r oc r and limited only by the amount of en- 
thalpy available at the sonic point, in our case to < 4. 

We have placed the location of the outflow boundaries 
out of causal contact with the propagating jet beam of 
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Figure 19. Mock observations based on model A at 43GHz and 86GHz for viewing angles i £ (0.01°, 20°, 40°) with resolution of 50/ias. 
Contours according to 43GHz intensity and spaced by factors of two. From left to right: Spectral index a, 43GHz polarization degree II 
and rotation measure in optically thin regions. 
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Figure 20. The dependence on viewing angle i of unresolved 
quantities for 86 GHz (red) and 43 GHz (blue). Top to bottom: 
Beamed flux, spectral index, integral polarization degree, observed 
polarization vector \i an d apparent rotation measure derived from 
Axf- 



interest. Thus, we can be sure that the calculated jet 
structure is purely self-collimated, and does not suffer 
from spurious boundary effects leading to an artificial 
collimation. We have investigated jets with a variety of 
poloidal electric current distributions. We find - some- 
what surprisingly - that the topology of the current dis- 
tribution, e.g. closed current circuits in comparison to 
current-carrying models, has little influence on the jet 
collimation. 

6.2. Jet radiation 

The benefit of having performed relativistic MHD sim- 
ulations of jet formation is that we could apply them to 
produce dynamically consistent emission maps to pre- 
dict VLBI radio and (sub-) mm observations of nearby 
AGN cores. For this purpose, we have developed a spe- 
cial relativistic synchrotron transport code fully taking 
into account self-absorption and internal Faraday rota- 
tion. Since the acceleration of non-thermal particles can 
not be followed self-consistently within the framework of 
pure MHD, it remains necessary to resort the particle en- 
ergy distribution to simple recipes. We have compared 
three prescriptions of the non-thermal particle energy 
distribution. We found good agreement in the alignment 
of the polarization structure, but considerable differences 



in the intensity maps. Thus, the polarization maps de- 
rived in this work can be considered as robust, while the 
intensities distribution should be regarded with caution. 

The str ict bi-modality of th e polarization direction sug- 
gested by Pariev et al. ( 2003 ) ; Lyutikov et al. ( 2005 ) and 
others can be circumvented when the structure ot a colli- 
mating jet is considered. However, the efficient collima- 
tion to near-cylindrical jet flows in general confirms these 
results obtained for optically thin cylindrical flows when 
the fast jet is considered. Thus, depending on the pitch 
angles of the emission region, also a spine-and-sheath po- 
larization structure could be observed. The relativistic 
swing effect skews the polarization compared to the non- 
re lativistic case. Our ra diati on models affirm the finding 
of Lyutikov et al.| ( 2005 ) and Clausen-Brown et al. ( 201 1 1 
that relativistic aberration promotes asymmetries m the 
polarization (half spine-sheaths) and also in the spec- 
tral index. The observational detection of such features 
would allow to determine the spin direction of the jet 
driver, be it the accretion disk or the central black hole. 

The frequency-dependent core shift in the radiation 
maps following our jet simulations is consisten t with an- 
alytical estimates of conical jets by |Lobanov ( 1998 1 in 



two jet models and slightly steeper in tfie third case con- 
sidered. We attribute this discrepancy to fact of jet colli- 
mation. The overall good agreement with the analytical 
estimate suggests that the standard diagnostics should 
provide robust results capable of determining the jet pa- 
rameters. With our radiation models we have confirmed 
the intuition that unresolved core shifts should lead to 
a breaking of the A 2 rotation measure law. Further, we 
have demonstrated that that law can be restored again 
as soon as the resolution is increased. Opacity effects do 
not allow to obtain a consistent A 2 -law across the spec- 
tral peak. Once the regimes are separated however, we 
obtain two valid relations for which the optically thin 
rotation measure is substantially increased over the op- 
tically thick case as it peers deeper into the Faraday pit. 

The interpretation of observations featuring both inter- 
nal Faraday rotation and changes in opacity is one of the 
most challenging aspects in polarimetric imaging of jets. 
With our detailed modeling, we were able to disentangle 
the depolarizing effects of opacity transition, differential 
Faraday rotation, and also beam effects such as ordinary 
beam depolarization or Faraday dispersion for two exem- 
plary jet models. We find that the unresolved, optically 
thin mm-wavelength radiation is depolarized due to both 
the changing emission geometry (down to ~ 40%), and 
the additional beam depolarization (down to ~ 30%). 
Increasing the Faraday opacity by observing at lower fre- 
quencies would lead to depolarization below the 1% level 
due to 1. Faraday dispersion and 2. differential Faraday 
rotation. 

We have also investigated the influence of resolution on 
the detectability of rotation measure (RM) gradients in 
the optically thin parsec-scale jet-core previewing mm- 
VLBI observations. To detect such gradients across the 
jet, a resolution of ~ 100rg would be required. Increas- 
ing the beam size leads to more monotonic transversal 
RM profiles. We find the peak magnitude of the RM to 
increase with resolution. High RMs beyond 10 4 radm -2 
are required to obtain a noticeable deflection in the mm- 
wavelength range. From the sources where high resolu- 
tion data is already available, namely Sgr A* and M87, 
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such high rotation measures are in fact observed, and we 
predict that many more objects in this class will be found 
at the advent of ALMA and global mm-VLBI. 

Finally, we have presented mock observations of spec- 
tral index, polarization degree and rotation measure for 
various inclinations. Asymmetries in the spectral index 
and polarization degree can be observed most clearly at 
high inclinations > 30°. The necessary resolution for 
this detection in a fiducial low Faraday rotation case (our 
model A) amounts to 50/ias, which could be reached with 
the next generation space- VLBI. At ~ 30°, the predom- 
inant polarization vector flips from perpendicular align- 
ment (with respect to the projected jet direction) for the 
blazar case to parallel alignment for the radio galaxy case 
at high inclinations. The flip in polarization is clearly 
detectable also from unresolved quantities. In summary, 
these findings suggest a bimodal distribution of the ob- 
served polarization direction of AGN core jets. However, 
by adding a substantial amount of Faraday rotation (our 
model B), this signature will be scrambled unless the ob- 
serving frequency is chosen high enough - confirming the 
popular intuition. 

In this paper we have focussed on general signatures of 
the synchrotron radiation in the large-scale helical fields 



in the acceleration region of relativistic MHD jets. With 
the developed tool set in hand, further progress can be 
made when calibrating the observational diagnostics with 
the mock observations that are detailed here. We expect 
a substantial improvement from a more consistent treat- 
ment of the non-thermal particles, taking into account 
particle acceleration and cooling. Potentially, also mod- 
eling of individual sources, or the cumulative statistics 
of AGN surveys applying dynamical simulations could 
be undertaken in the future. 
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APPENDIX 
STOKES TRANSPORT 



In the {I,I r , U } basis, the linear transfer equation dl/dl = £ — A I (20) has the coefficients 
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Naturally, e^ e ^ and K^ e,f> ' indicate the comoving emissivity/absorptivity in direction e' = n' x B' respectively 
b' = n' x e' given by line of sight n' and magnetic field direction B', while k — 1/2(k^ + k^)) is the average 
absorption coefficient. We apply the standa rd expressions for e^ e,b ^ and K^ e,b ^ valid for isotropic power-law particle 
distributions with index p — 2a + 1 following Pacholczyk ( 1970 ). For completeness, the expressions read: 
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depending on the constants 
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in Gauss cgs units where T(-) denotes the Gamma-function. The upper sign corresponds to the direction of the main 
polarization axis (§'). 

While the transformation of line of sight n leads to the well known relativistic aberration, 



n' = Dn (D + l) 



7 



7 + 1 







(A16) 



which supplies for the angle between co- moving magnetic field and line of sight B' sind' = | n' x B'|, the change 



betwe en e' and its observed counterpart e introduces a relativistic "swing" of the polarization (Blandford & Konigl 
|1979| . The observer system quantities follow as 



n x q 



Vq 2 - (n • q) : 



q = B + nx(/3xB) 



(A17) 



according to the short formulation discovered by Lyutikov et al. (2003). Hence the local observer-system angle of the 
emission Xe as measured from the image axis 1 (the projection of jet axis onto the plane of the sky) is given by 



COSXe 



le 



sin Xe = n • 



1 x e 



(A18) 



and we can conveniently take advantage of the identity sin 2x e = 2 sin Xe cos Xe to calcu late sin 2% e that appears in 
the absorption matrix. The description of Faraday rotation in the observers system, eq. (21) directly enters into the 
transfer coefficients ai3 and 023. 

With these considerations we solve the ordinary differential equation ( 20 1 using a fourth-order Runge-Kutta scheme 
for a grid of lines of sight. In practice this is realized by transforming the \r,z) plane information of the simulations into 
a coarsened cartesian grid with typically 100 x 100 x 500 cells by means of a Delaunay triangulation. For numerical 
stability, the adaptive step size is chosen small enough to satisfy Al < 6rg, At < 0.5 and Axf < tt/16 yielding 
sufficient convergence of the solution. 
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